#include <math.h>
#include <stdio.h>
#define  ROOT2       1.414213562

void get_ij_k_rands_(long long int *, long long int *, long long int *, 
		     float *, float *);
void amp_phase_(int *, int *, int *, double *, double *);
void gvar_(int *, int *, int *, float *, float *);
void code_mode_(int *,int *, int *, long long int *);
void principal_(int *, int *, int *, int *);
float ran4_(int *);
void self_test(void);
void setrand_(int *);



void get_ij_k_rands_(long long int *nft, long long int *ii, long long int *jj, 
		     float *ggvar1, float *ggvar2)
{
  int i,veclen;
  int ix,iy,iz;
  veclen = *nft/2;
  double ampl,phase;
  float g1,g2;
  int isign;

  ix = *ii;
  iy = *jj;
  if (ix>veclen) ix=ix-*nft;
  if (iy>veclen) iy=iy-*nft;
  for (i=0; i<veclen;i++){
    iz = i;
    if (iz>veclen) iz=iz-*nft;
    gvar_(&ix,&iy,&iz,&ggvar1[i],&ggvar2[i]);
  }	 
}

void amp_phase_(int *lq, int *mq, int *nq, double *ampl, double *phase)
{
const long long int iconst1 = 1003003002;
const long long int iconst2 = 1003003999;
int idmod,idmod2,idmod3,idmod4,irat1,irat2;
const float twopi = 2.*3.1415926535;
const float root2 = 1.414213562;
int iranA, iranB, iranC, iranD;
float ra1,ra2,ra3,ra4,amp_exp;
int l,m,n,ll,mm,nn,isign,iamp;
long long int id,id2;


 iamp = 1;  /* 1=Rayleigh dist   0=mean power */ 
 l = *lq;
 m = *mq;
 n = *nq;

 principal_(&l,&m,&n,&isign);

 if (isign==0){
      *ampl = 0.0;    /* DC mode */
      *phase= 0.0;
      return;
 }
 


 ll = -l;
 mm = -m;
 nn = -n;

 if (isign==1){
    code_mode_(&l, &m, &n, &id);
    code_mode_(&ll,&mm,&nn,&id2);
 }


 if (isign==-1){
    code_mode_(&l, &m, &n, &id2);
    code_mode_(&ll,&mm,&nn,&id);
 }

 idmod = (int)(id%iconst1);
 idmod2= (int)(id2%iconst1);
 irat1 = (int)((id/iconst1));
 irat2 = (int)((id2/iconst1));

 if ((irat1==0)&&(irat2==0)){


        iranA = 2*idmod-1;
        iranB = 2*idmod;
        iranC = 2*idmod2-1;
        iranD = 2*idmod2;


        ra1 = ran4_(&iranA);
        ra2 = 0.0;
        ra3 = ran4_(&iranC);
        ra4 = 0.0;
 }else{
   if ((irat1!=0)&&(irat2!=0)){

        idmod3 = (int)(id%iconst2);
        idmod4 = (int)(id2%iconst2);

        iranA = idmod;
        iranB = idmod2;
        iranC = idmod4;
        iranD = idmod3;
  

        ra1 = ran4_(&iranA);
        ra2 = ran4_(&iranB);
        ra3 = ran4_(&iranC);
        ra4 = ran4_(&iranD);

   }else{
           printf("Should not reach this point!\n");
           abort();
   }

 }

        ra1 = ra1 + ra4;
        if (ra1>1.0) ra1 = ra1-1.0;
        ra3 = ra3 + ra2;
        if (ra3>1.0) ra3 = ra3-1.0;
        if (ra3<=0.0) ra3 = ra3 + 1.e-8;
 

 
 if (iamp==0){
     *ampl = 1.0;
     *phase = -(double)(twopi*ra1);
 }else{

     *ampl = (double) ra3;
     *phase = -(double)(twopi*ra1);  
  }

      if ((l>0)&&(n<0)) *phase = -*phase;
      if ((l<0)&&(n>0)) *phase = -*phase;

      if ((l==0)&&(  ((m>0)&&(n<0)) || ((m<0)&&(n>0)))) *phase = -*phase;
      return;
}


void gvar_(int *lq, int *mq, int *nq, float *gvar1, float *gvar2)
{
const long long int iconst1 = 1003003002;
const long long int iconst2 = 1003003999;
int idmod,idmod2,idmod3,idmod4,irat1,irat2;
const float twopi = 2.*3.1415926535;
const float root2 = 1.414213562;
int iranA, iranB, iranC, iranD;
float ra1,ra2,ra3,ra4,amp_exp;
int l,m,n,ll,mm,nn,isign,iamp;
long long int id,id2;


 iamp = 1;  /* 1=Rayleigh dist   0=mean power */ 
 l = *lq;
 m = *mq;
 n = *nq;

 principal_(&l,&m,&n,&isign);

 if (isign==0){
      *gvar1 = 0.;    /* DC mode */
      *gvar2 = 0.;
      return;
 }
 


 ll = -l;
 mm = -m;
 nn = -n;

 if (isign==1){
    code_mode_(&l, &m, &n, &id);
    code_mode_(&ll,&mm,&nn,&id2);
 }


 if (isign==-1){
    code_mode_(&l, &m, &n, &id2);
    code_mode_(&ll,&mm,&nn,&id);
 }

 idmod = (int)(id%iconst1);
 idmod2= (int)(id2%iconst1);
 irat1 = (int)((id/iconst1));
 irat2 = (int)((id2/iconst1));

 if ((irat1==0)&&(irat2==0)){


        iranA = 2*idmod-1;
        iranB = 2*idmod;
        iranC = 2*idmod2-1;
        iranD = 2*idmod2;


        ra1 = ran4_(&iranA);
        ra2 = 0.0;
        ra3 = ran4_(&iranC);
        ra4 = 0.0;
 }else{
   if ((irat1!=0)&&(irat2!=0)){

        idmod3 = (int)(id%iconst2);
        idmod4 = (int)(id2%iconst2);

        iranA = idmod;
        iranB = idmod2;
        iranC = idmod4;
        iranD = idmod3;
      

        ra1 = ran4_(&iranA);
        ra2 = ran4_(&iranB);
        ra3 = ran4_(&iranC);
        ra4 = ran4_(&iranD);

   }else{
           printf("Should not reach this point!\n");
           abort();
   }

 }

        ra1 = ra1 + ra4;
        if (ra1>1.0) ra1 = ra1-1.0;
        ra3 = ra3 + ra2;
        if (ra3>1.0) ra3 = ra3-1.0;
        if (ra3<=0.0) ra3 = ra3 + 1.e-8;
 

 
 if (iamp==0){
     *gvar1 = cos(twopi*ra1)*root2;
     *gvar2 = (float)isign*sin(twopi*ra1)*root2;
      return;
 }else{
   
     *gvar1 = cos(twopi*ra1)*root2;
     *gvar2 = (float)isign * sin(twopi*ra1)*root2;

     
     amp_exp = sqrt(-log(ra3));


     *gvar1 = *gvar1 * amp_exp;
     *gvar2 = *gvar2 * amp_exp;
     return;  
      }

}





void setrand_(int *ns1)
{
 float rub;
 int idum,nrseed;

 self_test();


 nrseed = *ns1;

 idum = -abs(nrseed);

 rub = ran4_(&idum);   

 /*printf("Using random seed %d\n",*ns1);*/

}


/* Routine to check portability */
void self_test(void)
{
 float rub;
 int idum;
 int l,m,n;
 float gvar1,gvar2;

 if (sizeof(int)!=4){
   printf("Unexpected size for int type"); 
   abort();
     };
 if (sizeof(long long int)!=8){
   printf("Unexpected size for long long int type"); 
   abort();
     };

 /*  Check routines return standard values for a selection of l,m,n values */

 idum = -1379384218;
 rub = ran4_(&idum);   

 l = 23; 
 m = 13;
 n = -13;

 gvar_(&l,&m,&n,&gvar1,&gvar2);

 rub = fabs(gvar1 + 0.17146677);
 if (rub>1.e-8){
    printf("Error in self_test (23,13,-13) gvar1 %f \n",rub);
    abort();
 }
 rub = fabs(gvar2 + 0.22499602);
 if (rub>1.e-8){
    printf("Error in self_test (23,13,-13) gvar2 %f \n",rub);
    abort();
 }


 l = 1142;
 m = -2324;
 n = 9545;

 gvar_(&l,&m,&n,&gvar1,&gvar2);

 rub = fabs(gvar1 + 1.2492099);
 if (rub>3.e-8){
    printf("Error in self_test (1142,-2324,9545) gvar1 %14.9f \n",rub);
    abort();
 }
 rub = fabs(gvar2 + 0.055993884);
 if (rub>3.e-8){
    printf("Error in self_test (1142,-2324,9545) gvar2 %f \n",rub);
    abort();
 }

}




void code_mode_(int *lq, int *mq, int *nq, long long int *id)
{
 

 long long int l,m,n,iran,ish,idmin;
 long long int ll,mm,nn;
 
 l = *lq;
 m = *mq;
 n = *nq;

 
 iran = (abs(*lq)>abs(*mq))? abs(*lq):abs(*mq);
 iran = ((int)iran>abs(*nq))? (int)iran:abs(*nq);

 ish = 2*iran + 1;

 ll = l + iran;
 mm = m + iran;
 nn = n + iran;

 
 if (iran==0){
     *id = 1; 
     return;
 }

 idmin = 1 + (ish-2) * (ish-2) * (ish-2);

 
 if (ll==0){
     *id = idmin + ish*mm + nn;
     return;
 }

 if (ll==(ish-1)){
     *id = idmin + 5*ish*ish - 12*ish + 8 + ish*mm + nn;
     return;
 }

 idmin = idmin + ish*ish + 4*(ish-1)*(ll-1);


 if (mm>=nn){
     *id = idmin + mm + nn;
     return;
 }else{
     *id = idmin + 2*ish - 2 + mm + nn;
 }

}


void principal_(int *l, int *m, int *n, int *isign)
{


  if (*l!=0){
    if (*l>0){
              *isign = 1;
    }else{
              *isign= -1;
    }
    return;
  }

  if (*m!=0){
    if (*m>0){
              *isign = 1;
    }else{
              *isign = -1;
    }
    return;
  }

  if (*n!=0){
    if (*n>0){
              *isign =  1;
    }else{
              *isign = -1;
    }
    return;
  }

  *isign = 0;  /* DC mode */

}



#define NITER 4
 
void psdes(unsigned int* lword, unsigned int *irword)
 {
  unsigned int i,ia,ib,iswap,itmph=0,itmpl=0;
  static unsigned int c1[NITER]={
  0xbaa96887L, 0x1e17d32cL, 0x03bcdc3cL, 0x0f33d1b2L};
  static unsigned int c2[NITER]={
  0x4b0f3b58L, 0xe874f0c3L, 0x6955c5a6L, 0x55a7ca46L};
             	
  for(i=0;i<NITER;i++)
    {
         ia=(iswap=(*irword)) ^ c1[i];
         itmpl = ia & 0xffff;
         itmph = ia >> 16;
         ib=itmpl*itmpl+ ~(itmph*itmph);
        *irword=*lword ^ (((ia=(ib>>16) | ((ib & 0xffff) << 16)) ^ c2[i])+itmpl*itmph);
	*lword=iswap;
           	}
             }

            
  float ran4_(int *idum)
   {
            	unsigned int irword,itemp,lword;
            	static int idums =0;
             

            	static unsigned int jflone = 0x3f800000;
            	static unsigned int jflmsk = 0x007fffff;
            
            	if(*idum<0)
           	{
            		idums = -(*idum);
            		*idum=1;
             	}
             
            	irword=(*idum);
            	lword=idums;
             	psdes(&lword,&irword);
             	itemp=jflone | (jflmsk&irword);
             	++(*idum);
            	return (*(float *)&itemp)-1.0;
             }
             
